GEOG-533 Lab 9

Question 1

Use R package UScensus2010county to complete the following tasks: (20 pt.)

Question 1(a)

Plot a map of New York counties using the plot function.

library(rgdal)
if (!require(UScensus2010county)) {install.county("osx")}
if (!require(UScensus2010tract)) {install.tract("osx")}
library(UScensus2010)
library(UScensus2010county)
data(new_york.county10)
shp <- new_york.county10
plot(shp)

Question 1(b)

Plot a map of New York counties using the qtm function.

library(tmap)
qtm(shp)

Question 1(c)

How many counties in New York State?

nrow(shp)

[1] 62

## There are 62 counties in New York State

Question 1(d)

What’s the 3-digit fips code of Broome County?

df <- shp@data
df

## The FIPS code of Broome County is 007

Question 1(e)

Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.

sum(df$P0010001)

[1] 19378102

min(df$P0010001)

[1] 4836

max(df$P0010001)

[1] 2504700

mean(df$P0010001)

[1] 312550

median(df$P0010001)

[1] 91301

library(moments)
skewness(df$P0010001)

[1] 2.579801

## There is a total population of 19378102. The minimum population of all the counties is 4836, the maximum is 2504700, the mean is 312550, the median is 91301, and the skewness is 2.579801.

Question 1(f)

Create a histogram and a boxplot of the population.

Population <- df$P0010001
hist(Population)

boxplot(Population)

Question 2

Use R package UScensus2010tract to complete the following tasks: (20 pt.)

library(UScensus2010tract)

Question 2(a)

Plot a map of New York census tracts using the plot function.

data(new_york.tract10)
shp <- new_york.tract10
plot(shp)

Question 2(b)

Compute the total population based on census tracts.

df <- shp@data
df

sum(df$P0010001)

[1] 19378102

## The total population based on census tracts is 19378102

Question 2(c)

Select all census tracts in Broome County and plot the map.

sel <- df$county == "007"
plot(shp[sel,])

Question 2(d)

What’s the total population of Broome County?

sum(df$P0010001[df$county == "007"])

[1] 200600

## The total population of Broome County is 200600

Question 2(e)

Create a histogram and a boxplot of population based on census tracts of Broome County.

broome_county_population <- (df$P0010001[df$county == "007"])
hist(broome_county_population)

boxplot(broome_county_population)

Question 2(f)

Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.

broome_county_shp <- subset(shp, county == "007")
broome_county_tracts <- broome_county_shp[, 1:5]
broome_county_tracts$PopulationRatio <- (broome_county_tracts$P0010001 / 200600)

Question 3

Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)

library(raster)
library(ncdf4)

Question 3(a)

Load the multi-band raster dataset “NDVI.nc” into R.

setwd("~/Desktop/533/Blank Labs/Lab_09")

The working directory was changed to /Users/deborahgabay/Desktop/533/Blank Labs/Lab_09 inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

ndvi.rb <- brick("NDVI.nc")

Question 3(b)

Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.

ndvi.rb

class       : RasterBrick 
dimensions  : 360, 720, 259200, 36  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/deborahgabay/Desktop/533/Blank Labs/Lab_09/NDVI.nc 
names       : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ... 
Date        : 2000-01-01, 2002-12-01 (min, max)
varname     : NDVI 

Question 3(c)

Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.

aggregate(ndvi.rb)

class       : RasterBrick 
dimensions  : 180, 360, 64800, 36  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ... 
min values  :           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0,           0, ... 
max values  :   0.9089250,   0.8816000,   0.8811250,   0.9275750,   0.9325000,   0.9033999,   0.9327000,   0.9038250,   0.8804000,   0.8941250,   0.9066000,   0.9554500,   0.9094500,   0.9066250,   0.8673500, ... 

mean(ndvi.rb)

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /private/var/folders/cq/rlf0yv_17254sx86cgg4yrz00000gn/T/RtmpgYpMs5/raster/r_tmp_2019-04-30_231415_5280_98713.grd 
names       : layer 
values      : 0, 0.8537528  (min, max)

max(ndvi.rb)

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.9681  (min, max)

Question 3(d)

Plot the maps of monthly NDVI of the year 2001

ndvi2001 <- ndvi.rb[[13:24]]
plot(ndvi2001)

Question 3(e)

Create histograms of monthly NDVI of the year 2001.

hist(ndvi2001, 1)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 2)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 3)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 4)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 5)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 6)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 7)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 8)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 9)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 10)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 11)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

hist(ndvi2001, 12)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

Question 3(f)

Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.

hist(ndvi2001, 7)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

values <- ndvi2001[50,50]
plot(as.vector(values), type = "b")